I want to evaluate the function \[ f(x, y) = \frac{1}{|0.5 - x^4 - y^4| + 0.1} \] which is not particularly well-behaved in that it has a discontinuity due to the absolute value, and it does not disappear at the boundaries. Evaluating the function through classical sparse grids is insufficient because ot requires too many points to capture the ridge. Adaptive sparse grids are better suited for this problem. This write-up demonstrates the potential of adaptive sparse grids.

Brute Force

First, I directly evaluate the function over the unit square, \([0, 1]^2\) using \(100^2\) points.

library('plotly')

data.classic <- read.csv('/home/ahyan/Dropbox/Housing Market and Leverage Cycle/Code/src/Sparse/CSG.csv', header = FALSE)

fig1 <- plot_ly(data.classic, x = ~V1, y = ~V2, z = ~V3, marker = list(color = ~V3, colorscale = 'Portland', showscale = TRUE, size = 3))

fig1 <- fig1 %>% layout(title = 'Direct Evaluation')

fig1 <- fig1 %>% layout(scene = list(xaxis = list(title = 'x'),
                                     yaxis = list(title = 'y'),
                                     zaxis = list(title = 'f(x, y)')))

fig1

Classical Sparse Grids

Now I evaluate the function on sparse grids of level 14. In two dimensions, that means almost 280,000 grid points.

fig2 <- plot_ly(data.classic, x = ~V1, y = ~V2, z = ~V4, marker = list(color = ~V4, colorscale = 'Portland', showscale = TRUE, size = 3))

fig2 <- fig2 %>% layout(title = 'Classical Sparse Grid')

fig2 <- fig2 %>% layout(scene = list(xaxis = list(title = 'x'),
                                     yaxis = list(title = 'y'),
                                     zaxis = list(title = 'f(x, y)')))

fig2

The error between brute force and sparse grids, calculated as the \(L_2\) norm is 0.0889.

\[ L_2 = \frac{1}{N} \bigg(\sum_{i = 1}^N |f(x_i) - u(x_i)|^2\bigg)^\frac{1}{2} \]

Adaptive Sparse Grids

Finally, I evaluate the function with adaptive sparse grids with level 5. In two dimensions, the algo starts with 257 bounds (boundary included) and iterates until the error is less than 0.1 which happens with 6571 grid points and the error is 0.0897.

data.adaptive <- read.csv('/home/ahyan/Dropbox/Housing Market and Leverage Cycle/Code/src/Sparse/ASG.csv', header = FALSE)

fig3 <- plot_ly(data.adaptive, x = ~V1, y = ~V2, z = ~V4, marker = list(color = ~V4, colorscale = 'Portland', showscale = TRUE, size = 3))

fig3 <- fig3 %>% layout(title = 'Adaptive Sparse Grids')

fig3 <- fig3 %>% layout(scene = list(xaxis = list(title = 'x'),
                                     yaxis = list(title = 'y'),
                                     zaxis = list(title = 'f(x, y)')))

fig3

Life Cycle Problem

Basic

data.life <- read.csv('/home/ahyan/Dropbox/Housing Market and Leverage Cycle/Code/src/Sparse/Life.csv', header = FALSE)

fig4 <- plot_ly(data.life, x = ~V1, y = ~V2, z = ~V3, marker = list(color = ~V3, colorscale = 'Portland', showscale = TRUE, size = 3))

fig4 <- fig4 %>% layout(title = 'Direct')

fig4 <- fig4 %>% layout(scene = list(xaxis = list(title = 'e'),
                                     yaxis = list(title = 'x'),
                                     zaxis = list(title = 'V(e, x)')))

fig4
fig5 <- plot_ly(data.life, x = ~V1, y = ~V2, z = ~V4, marker = list(color = ~V4, colorscale = 'Portland', showscale = TRUE, size = 3))

fig5 <- fig5 %>% layout(title = 'Adaptive')

fig5 <- fig5 %>% layout(scene = list(xaxis = list(title = 'e'),
                                     yaxis = list(title = 'x'),
                                     zaxis = list(title = 'V(e, x)')))

fig5

With Nelder-Mead

data.lifeB <- read.csv('/home/ahyan/Dropbox/Housing Market and Leverage Cycle/Code/src/Sparse/LifeB.csv', header = FALSE)

fig6 <- plot_ly(data.lifeB, x = ~V1, y = ~V2, z = ~V3, marker = list(color = ~V3, colorscale = 'Portland', showscale = TRUE, size = 3))

fig6 <- fig6 %>% layout(title = 'Direct')

fig6 <- fig6 %>% layout(scene = list(xaxis = list(title = 'e'),
                                     yaxis = list(title = 'x'),
                                     zaxis = list(title = 'V(e, x)')))

fig6
LS0tCnRpdGxlOiAiU3BhcnNlIEdyaWRzIgphdXRob3I6IEFoeWFuIFBhbmp3YW5pCmRhdGU6IEF1Z3VzdCAxMXRoLCAyMDIwCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KCkkgd2FudCB0byBldmFsdWF0ZSB0aGUgZnVuY3Rpb24KJCQKZih4LCB5KSA9IFxmcmFjezF9e3wwLjUgLSB4XjQgLSB5XjR8ICsgMC4xfQokJAp3aGljaCBpcyBub3QgcGFydGljdWxhcmx5IF93ZWxsLWJlaGF2ZWRfIGluIHRoYXQgaXQgaGFzIGEgZGlzY29udGludWl0eSBkdWUgdG8gdGhlIGFic29sdXRlIHZhbHVlLCBhbmQKaXQgZG9lcyBub3QgZGlzYXBwZWFyIGF0IHRoZSBib3VuZGFyaWVzLiBFdmFsdWF0aW5nIHRoZSBmdW5jdGlvbiB0aHJvdWdoIGNsYXNzaWNhbCBzcGFyc2UgZ3JpZHMgaXMgaW5zdWZmaWNpZW50CmJlY2F1c2Ugb3QgcmVxdWlyZXMgdG9vIG1hbnkgcG9pbnRzIHRvIGNhcHR1cmUgdGhlIHJpZGdlLiBBZGFwdGl2ZSBzcGFyc2UgZ3JpZHMgYXJlIGJldHRlciBzdWl0ZWQgZm9yIHRoaXMKcHJvYmxlbS4gVGhpcyB3cml0ZS11cCBkZW1vbnN0cmF0ZXMgdGhlIHBvdGVudGlhbCBvZiBhZGFwdGl2ZSBzcGFyc2UgZ3JpZHMuCgojIyBCcnV0ZSBGb3JjZQoKRmlyc3QsIEkgZGlyZWN0bHkgZXZhbHVhdGUgdGhlIGZ1bmN0aW9uIG92ZXIgdGhlIHVuaXQgc3F1YXJlLCAkWzAsIDFdXjIkIHVzaW5nICQxMDBeMiQgcG9pbnRzLgoKYGBge3Igd2FybmluZz1GQUxTRSwgbWVzc2FnZT1GQUxTRSwgZmlnLmFsaWduPSdjZW50ZXInfQoKbGlicmFyeSgncGxvdGx5JykKCmRhdGEuY2xhc3NpYyA8LSByZWFkLmNzdignL2hvbWUvYWh5YW4vRHJvcGJveC9Ib3VzaW5nIE1hcmtldCBhbmQgTGV2ZXJhZ2UgQ3ljbGUvQ29kZS9zcmMvU3BhcnNlL0NTRy5jc3YnLCBoZWFkZXIgPSBGQUxTRSkKCmZpZzEgPC0gcGxvdF9seShkYXRhLmNsYXNzaWMsIHggPSB+VjEsIHkgPSB+VjIsIHogPSB+VjMsIG1hcmtlciA9IGxpc3QoY29sb3IgPSB+VjMsIGNvbG9yc2NhbGUgPSAnUG9ydGxhbmQnLCBzaG93c2NhbGUgPSBUUlVFLCBzaXplID0gMykpCgpmaWcxIDwtIGZpZzEgJT4lIGxheW91dCh0aXRsZSA9ICdEaXJlY3QgRXZhbHVhdGlvbicpCgpmaWcxIDwtIGZpZzEgJT4lIGxheW91dChzY2VuZSA9IGxpc3QoeGF4aXMgPSBsaXN0KHRpdGxlID0gJ3gnKSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHlheGlzID0gbGlzdCh0aXRsZSA9ICd5JyksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB6YXhpcyA9IGxpc3QodGl0bGUgPSAnZih4LCB5KScpKSkKCmZpZzEKCmBgYAoKIyMgQ2xhc3NpY2FsIFNwYXJzZSBHcmlkcwoKTm93IEkgZXZhbHVhdGUgdGhlIGZ1bmN0aW9uIG9uIHNwYXJzZSBncmlkcyBvZiBsZXZlbCAxNC4gSW4gdHdvIGRpbWVuc2lvbnMsIHRoYXQgbWVhbnMgYWxtb3N0IDI4MCwwMDAgZ3JpZCBwb2ludHMuCgpgYGB7ciB3YXJuaW5nPUZBTFNFLCBtZXNzYWdlPUZBTFNFLCBmaWcuYWxpZ249J2NlbnRlcid9CmZpZzIgPC0gcGxvdF9seShkYXRhLmNsYXNzaWMsIHggPSB+VjEsIHkgPSB+VjIsIHogPSB+VjQsIG1hcmtlciA9IGxpc3QoY29sb3IgPSB+VjQsIGNvbG9yc2NhbGUgPSAnUG9ydGxhbmQnLCBzaG93c2NhbGUgPSBUUlVFLCBzaXplID0gMykpCgpmaWcyIDwtIGZpZzIgJT4lIGxheW91dCh0aXRsZSA9ICdDbGFzc2ljYWwgU3BhcnNlIEdyaWQnKQoKZmlnMiA8LSBmaWcyICU+JSBsYXlvdXQoc2NlbmUgPSBsaXN0KHhheGlzID0gbGlzdCh0aXRsZSA9ICd4JyksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB5YXhpcyA9IGxpc3QodGl0bGUgPSAneScpLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgemF4aXMgPSBsaXN0KHRpdGxlID0gJ2YoeCwgeSknKSkpCgpmaWcyCmBgYAoKVGhlIGVycm9yIGJldHdlZW4gYnJ1dGUgZm9yY2UgYW5kIHNwYXJzZSBncmlkcywgY2FsY3VsYXRlZCBhcyB0aGUgJExfMiQgbm9ybSBpcyAwLjA4ODkuCgokJApMXzIgPSBcZnJhY3sxfXtOfSBcYmlnZyhcc3VtX3tpID0gMX1eTiB8Zih4X2kpIC0gdSh4X2kpfF4yXGJpZ2cpXlxmcmFjezF9ezJ9CiQkCgojIyBBZGFwdGl2ZSBTcGFyc2UgR3JpZHMKCkZpbmFsbHksIEkgZXZhbHVhdGUgdGhlIGZ1bmN0aW9uIHdpdGggYWRhcHRpdmUgc3BhcnNlIGdyaWRzIHdpdGggbGV2ZWwgNS4gSW4gdHdvIGRpbWVuc2lvbnMsIHRoZSBhbGdvIHN0YXJ0cwp3aXRoIDI1NyBib3VuZHMgKGJvdW5kYXJ5IGluY2x1ZGVkKSBhbmQgaXRlcmF0ZXMgdW50aWwgdGhlIGVycm9yIGlzIGxlc3MgdGhhbiAwLjEgd2hpY2ggaGFwcGVucyB3aXRoIDY1NzEKZ3JpZCBwb2ludHMgYW5kIHRoZSBlcnJvciBpcyAwLjA4OTcuCgpgYGB7ciB3YXJuaW5nPUZBTFNFLCBtZXNzYWdlPUZBTFNFLCBmaWcuYWxpZ249J2NlbnRlcid9CmRhdGEuYWRhcHRpdmUgPC0gcmVhZC5jc3YoJy9ob21lL2FoeWFuL0Ryb3Bib3gvSG91c2luZyBNYXJrZXQgYW5kIExldmVyYWdlIEN5Y2xlL0NvZGUvc3JjL1NwYXJzZS9BU0cuY3N2JywgaGVhZGVyID0gRkFMU0UpCgpmaWczIDwtIHBsb3RfbHkoZGF0YS5hZGFwdGl2ZSwgeCA9IH5WMSwgeSA9IH5WMiwgeiA9IH5WNCwgbWFya2VyID0gbGlzdChjb2xvciA9IH5WNCwgY29sb3JzY2FsZSA9ICdQb3J0bGFuZCcsIHNob3dzY2FsZSA9IFRSVUUsIHNpemUgPSAzKSkKCmZpZzMgPC0gZmlnMyAlPiUgbGF5b3V0KHRpdGxlID0gJ0FkYXB0aXZlIFNwYXJzZSBHcmlkcycpCgpmaWczIDwtIGZpZzMgJT4lIGxheW91dChzY2VuZSA9IGxpc3QoeGF4aXMgPSBsaXN0KHRpdGxlID0gJ3gnKSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHlheGlzID0gbGlzdCh0aXRsZSA9ICd5JyksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB6YXhpcyA9IGxpc3QodGl0bGUgPSAnZih4LCB5KScpKSkKCmZpZzMKYGBgCgoKIyMgTGlmZSBDeWNsZSBQcm9ibGVtCgpCYXNpYwpgYGB7ciB3YXJuaW5nPUZBTFNFLCBtZXNzYWdlPUZBTFNFLCBmaWcuYWxpZ249J2NlbnRlcid9CgpkYXRhLmxpZmUgPC0gcmVhZC5jc3YoJy9ob21lL2FoeWFuL0Ryb3Bib3gvSG91c2luZyBNYXJrZXQgYW5kIExldmVyYWdlIEN5Y2xlL0NvZGUvc3JjL1NwYXJzZS9MaWZlLmNzdicsIGhlYWRlciA9IEZBTFNFKQoKZmlnNCA8LSBwbG90X2x5KGRhdGEubGlmZSwgeCA9IH5WMSwgeSA9IH5WMiwgeiA9IH5WMywgbWFya2VyID0gbGlzdChjb2xvciA9IH5WMywgY29sb3JzY2FsZSA9ICdQb3J0bGFuZCcsIHNob3dzY2FsZSA9IFRSVUUsIHNpemUgPSAzKSkKCmZpZzQgPC0gZmlnNCAlPiUgbGF5b3V0KHRpdGxlID0gJ0RpcmVjdCcpCgpmaWc0IDwtIGZpZzQgJT4lIGxheW91dChzY2VuZSA9IGxpc3QoeGF4aXMgPSBsaXN0KHRpdGxlID0gJ2UnKSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHlheGlzID0gbGlzdCh0aXRsZSA9ICd4JyksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB6YXhpcyA9IGxpc3QodGl0bGUgPSAnVihlLCB4KScpKSkKCmZpZzQKCmBgYAoKCgpgYGB7ciB3YXJuaW5nPUZBTFNFLCBtZXNzYWdlPUZBTFNFLCBmaWcuYWxpZ249J2NlbnRlcid9CgpmaWc1IDwtIHBsb3RfbHkoZGF0YS5saWZlLCB4ID0gflYxLCB5ID0gflYyLCB6ID0gflY0LCBtYXJrZXIgPSBsaXN0KGNvbG9yID0gflY0LCBjb2xvcnNjYWxlID0gJ1BvcnRsYW5kJywgc2hvd3NjYWxlID0gVFJVRSwgc2l6ZSA9IDMpKQoKZmlnNSA8LSBmaWc1ICU+JSBsYXlvdXQodGl0bGUgPSAnQWRhcHRpdmUnKQoKZmlnNSA8LSBmaWc1ICU+JSBsYXlvdXQoc2NlbmUgPSBsaXN0KHhheGlzID0gbGlzdCh0aXRsZSA9ICdlJyksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB5YXhpcyA9IGxpc3QodGl0bGUgPSAneCcpLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgemF4aXMgPSBsaXN0KHRpdGxlID0gJ1YoZSwgeCknKSkpCgpmaWc1CmBgYAoKV2l0aCBOZWxkZXItTWVhZApgYGB7ciB3YXJuaW5nPUZBTFNFLCBtZXNzYWdlPUZBTFNFLCBmaWcuYWxpZ249J2NlbnRlcid9CgpkYXRhLmxpZmVCIDwtIHJlYWQuY3N2KCcvaG9tZS9haHlhbi9Ecm9wYm94L0hvdXNpbmcgTWFya2V0IGFuZCBMZXZlcmFnZSBDeWNsZS9Db2RlL3NyYy9TcGFyc2UvTGlmZUIuY3N2JywgaGVhZGVyID0gRkFMU0UpCgpmaWc2IDwtIHBsb3RfbHkoZGF0YS5saWZlQiwgeCA9IH5WMSwgeSA9IH5WMiwgeiA9IH5WMywgbWFya2VyID0gbGlzdChjb2xvciA9IH5WMywgY29sb3JzY2FsZSA9ICdQb3J0bGFuZCcsIHNob3dzY2FsZSA9IFRSVUUsIHNpemUgPSAzKSkKCmZpZzYgPC0gZmlnNiAlPiUgbGF5b3V0KHRpdGxlID0gJ0RpcmVjdCcpCgpmaWc2IDwtIGZpZzYgJT4lIGxheW91dChzY2VuZSA9IGxpc3QoeGF4aXMgPSBsaXN0KHRpdGxlID0gJ2UnKSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHlheGlzID0gbGlzdCh0aXRsZSA9ICd4JyksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB6YXhpcyA9IGxpc3QodGl0bGUgPSAnVihlLCB4KScpKSkKCmZpZzYKCmBgYAo=